% Script calculates the solution of integral equations for the current
% distribution and input impedance of the antenna of radius a located on
% the surface of plasma column of radius b in the case of resonant
% plasma inside the column

clear all;
clc;
close all;
format short e;

% set constants
CONSTS = [];
% f
% (18e+6 Hz --- 400e+6 Hz @ B_0 = 6)
% (30e+6 Hz --- 1200e+6 Hz @ B_0 = 10)
CONSTS.f = 27.12e+6; %0.1e+9; % 
CONSTS.c = 3e+10; % cm/s
CONSTS.a = 0.5; % cm
CONSTS.gamma = 0.577;
CONSTS.e = -4.80320427e-10;
CONSTS.m_e = 9.10660e-28;
CONSTS.B_0 = 800; %200; % 620 200 5 --- 800 gauss 200
CONSTS.N = 1e+13; %1e+11; % 1e+11 --- 1e+13 cm^-3
CONSTS.d = (CONSTS.a)/50; % 0.1 cm --- 0,01 mm
CONSTS.delta = 0.05; % 4.39e-2 rad for a = 0.4 --- 1.9 cm
CONSTS.phi_0 = 0;
CONSTS.eta_a = 1;
CONSTS.use_S = true;

CONSTS = update_CONSTS(CONSTS);

% if(true)
%     color_vec = ['b', 'r', 'g'];
%     figure;
%     hold;
%     title('abs I_{\Sigma}(x)');
%     xlabel('x, cm');
%     ylabel('(|I_{\Sigma}(x)/I_{\Sigma max}|)');
%     legend('(I(x))', ...
%         'location', 'SouthEast');
%     h_abs = gca; 
% 
%     figure;
%     hold;
%     title('real part I_{\Sigma}(x)');
%     xlabel('x, cm');
%     ylabel('Re(I_{\Sigma}(x)/I_{\Sigma max})');
%     legend('Re(I(x))', ...
%         'location', 'SouthEast');
%     h_re = gca; 
% 
%     figure;
%     hold;
%     title('imag part I_{\Sigma}(x)');
%     xlabel('x, cm');
%     ylabel('Im(I_{\Sigma}(x)/I_{\Sigma max})');
%     legend('Im(I(x))', ...
%         'location', 'SouthEast');
%     h_im = gca; 
% 
%     figure;
%     hold;
%     title('phase I_{\Sigma}(x)');
%     xlabel('x, cm');
%     ylabel('angle(I_{\Sigma}(x))');
%     legend('angle(I(x))', ...
%         'location', 'SouthEast');
%     h_phase = gca; 
% end
% 
% c = CONSTS.c;
% Z0=4*pi/c;
% k0 = CONSTS.k0;
% d = CONSTS.d;
% coef_sgs_to_si = 9e+11;
% 
% x = (0:0.1:50)';
% 
% N_vec = linspace(1e+11, 1e+13, 3)';
% % y = zeros(size(N_vec),size(x));
% for i=1:size(N_vec,1)
%     N=N_vec(i);
%     LOCAL_CONSTS = CONSTS;
%     LOCAL_CONSTS.N = N;
%     LOCAL_CONSTS = update_CONSTS(LOCAL_CONSTS);
%     eps = LOCAL_CONSTS.eps;
%     eta = LOCAL_CONSTS.eta;
% 
%     fprintf('eps = %e\n', eps);
%     fprintf('eta = %e\n', eta);
%     
%     eps_eff = ((1+abs(eps*eta))*sqrt(abs(eps*eta)))/(2*(1+abs(eps*eta)+1i*sqrt(abs(eps*eta))));
%     h = -k0*sqrt(-1i*eps_eff);
%     h_odn = k0*(abs(eps*eta))^(1/4)*(1-1i)/sqrt(2);
%     h_fs = k0;
%     y = -((pi*h)/(Z0*k0*log(4/(k0*d))))*exp(-1i*h*abs(x));
%     y_odn = (pi*h_odn/Z0*k0*log(4/(k0*d)))*exp(-1i*h_odn*abs(x));
%     y_fs = (pi*h_fs/Z0*k0*log(4/(k0*d)))*exp(-1i*h_fs*abs(x));
%     z(i) = -coef_sgs_to_si*(1i*Z0*k0/(pi*h))*log(4/(k0*d));
%     ymaxre = max(real(y));
%     ymaxim = max(imag(y));
%     ymaxabs = max(abs(y));
%     if(true)
%         plot(h_abs, x, abs(y./y(1)), [color_vec(i), '-'], ...
%                     x, abs(y_odn./y_odn(1)), [color_vec(i), '--'], ...
%                     x, abs(y_fs./y_fs(1)), [color_vec(i), 'o-']);
%         plot(h_re, x, real(y)./real(y(1)), [color_vec(i), '-'], x, real(y_odn)./real(y_odn(1)), [color_vec(i), '--']);
%         plot(h_im, x, imag(y)./imag(y(1)), [color_vec(i), '-'], x, imag(y_odn)./imag(y_odn(1)), [color_vec(i), '--']);
%         plot(h_phase, x, to_degrees(angle(y)), [color_vec(i), '-']);
%     end
% end

% figure; plot(x, abs(y./ymax), 'b');
% figure; plot(x, real(y./ymax), 'b');
% figure; plot(x, imag(y./ymax), 'b');
% figure; plot(N_vec, real(z), 'b', N_vec, imag(z), 'r-');


% P = CONSTS.P;
% P_c = CONSTS.P_c;
% plot_data = false;
% nx = 10;
% nz_vec = (1:0.1:2*P)';
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %------------------ nxyk --------------------------------------------------
% [nxy1 nxy2] = func_nxyk_from_nz(nz_vec, CONSTS); 
% if(plot_data)
%     figure; plot(nz_vec, real(nxy1), 'b-', nz_vec, real(nxy2), 'r.-' ...
%                  , [P_c, P_c], [0, max(real(nxy1))], 'm:' ...
%                  , [P, P], [0, max(real(nxy1))], 'k:' ...
%                  );
%     figure; plot(nz_vec, imag(nxy1), 'b-', nz_vec, imag(nxy2), 'r.-' ...
%                  , [P_c, P_c], [0, max(real(nxy1))], 'm:' ...
%                  , [P, P], [0, max(real(nxy1))], 'k:' ...
%                  );
%     figure; plot(nz_vec, abs(nxy1), 'b-', nz_vec, abs(nxy2), 'r.-');
% end
% 
% nxy = func_nxy_from_nz(nz_vec, CONSTS);
% if(plot_data)
%     figure; plot(nz_vec, real(nxy), 'b-', nz_vec, imag(nxy), 'r-' ...
%                  , [P_c, P_c], [0, max(abs(nxy))], 'm:' ...
%                  ... %, [P, P], [0, max(real(nxy1))], 'k:' ...
%                  );
% end
% 
% %------------------ delta -------------------------------------------------
% delta = func_delta_from_nz(nz_vec, nx, CONSTS);
% if(plot_data)
%     figure; plot(nz_vec, real(delta), 'b-', nz_vec, imag(delta), 'r-' ...
%                   ,[P_c, P_c], [0, max(abs(delta))], 'm:' ...
%                  ,[P, P], [0, max(abs(delta))], 'k:' ...
%                  );
%     figure; plot(nz_vec, abs(delta), 'b-' ...
%                  ,[P_c, P_c], [0, max(abs(delta))], 'm:' ...
%                  ,[P, P], [0, max(abs(delta))], 'k:' ...
%                  );
% end

%------------------ zeros of delta ----------------------------------------
% span = [1 P];
% [p_m coef] = func_poles(nx, CONSTS, span, plot_data);

nx = 10;
d = CONSTS.d;
nzstep = 0.05;
nzhigh = 1000;
zeta_step = 0.01;
nz = [(-nzhigh:nzstep:-1-nzstep) (-1+nzstep:nzstep:0-nzstep) (0+nzstep:nzstep:1-nzstep) (1+nzstep:nzstep:nzhigh)]';
zeta_vec = (zeta_step:zeta_step:1)'.*d;
if(false)
for i=1:size(zeta_vec,1)
    zeta = zeta_vec(i);
    if(true)
        re1 = quad_gauss(@(nz)(real(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [-nzhigh, -1-1/nzhigh], 0.005, 6);
        re2 = quad_gauss(@(nz)(real(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [-1+1/nzhigh, 0-1/nzhigh], 0.005, 6);
        re3 = quad_gauss(@(nz)(real(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [0+1/nzhigh, 1-1/nzhigh], 0.005, 6);
        re4 = quad_gauss(@(nz)(real(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [1+1/nzhigh, nzhigh], 0.005, 6);
        im1 = quad_gauss(@(nz)(imag(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [-nzhigh, -1-1/nzhigh], 0.005, 6);
        im2 = quad_gauss(@(nz)(imag(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [-1+1/nzhigh, 0-1/nzhigh], 0.005, 6);
        im3 = quad_gauss(@(nz)(imag(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [0+1/nzhigh, 1-1/nzhigh], 0.005, 6);
        im4 = quad_gauss(@(nz)(imag(funcKxfull(nx, nz, zeta, CONSTS))), ...
                            [1+1/nzhigh, nzhigh], 0.005, 6);
        Kxre(i) = re1+re2+re3+re4;
        Kxim(i) = (im1+im2+im3+im4);
    end
    if(false)
        Kxre(i) = trapz(nz, (real(funcKxfull(nx, nz, zeta, CONSTS))));
        Kxim(i) = trapz(nz, (imag(funcKxfull(nx, nz, zeta, CONSTS))));
        Kx1(i) = Kxre(i)+1i*Kxim(i);
    end


end
end
%%    
y = funcKxfull(nx, nz, 0.5, CONSTS);

figure; plot(nz, real(y), 'b', nz, imag(y), 'r');
% figure; plot(zeta_vec./d, Kxre, 'b', zeta_vec./d, Kxim, 'r');







